arXiv:1503.08926v2 [math.NA] 3 Sep 2015 


Convergence of summation-by-parts finite 
difference methods for the wave equation 


Siyang Wang*'!’ Gunilla Kreiss^ 


Abstract 

In this paper, we consider finite difference approximations of the sec¬ 
ond order wave equation. We use finite difference operators satisfying 
the summation-by-parts property to discretize the equation in space. 
Boundary conditions and grid interface conditions are imposed by the 
simultaneous-approximation-term technique. Typically, the truncation 
error is larger at the grid points near a boundary or grid interface than 
that in the interior. Normal mode analysis can be used to analyze how the 
large truncation error affects the convergence rate of the underlying stable 
numerical scheme. If the semi-discretized equation satisfies a determinant 
condition, two orders are gained from the large truncation error. However, 
many interesting second order equations do not satisfy the determinant 
condition. We then carefully analyze the solution of the boundary system 
to derive a sharp estimate for the error in the solution and acquire the 
gain in convergence rate. The result shows that stability does not auto¬ 
matically yield a gain of two orders in convergence rate. The accuracy 
analysis is verified by numerical experiments. 

Keywords: Second order wave equation, SBP-SAT, Finite difference, Ac¬ 
curacy, Convergence, Determinant condition, Normal mode analysis 
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1 Introduction 

In many physical problems, such as acoustics, seismology, and electromagnetism, 
the underlying equations can be formulated as systems of second order time de¬ 
pendent hyperbolic partial differential equations. For wave propagation prob¬ 
lems, it has been shown that high order accurate time marching methods as 
well as high order spatial discretization are more efficient to solve these prob¬ 
lems on smooth domains [ 319 ]. The major difficulties with high order spatial 
discretization are the numerical treatment of boundary conditions, grid interface 
conditions and interface conditions at material interfaces. 

High order finite difference methods have been widely used for wave propa¬ 
gation problems. We use summation-by-parts (SBP) finite difference operators 
[ni m na to approximate spatial derivatives. In order to guarantee strict 
stability and maintain high accuracy, boundary conditions and grid interface 
conditions are imposed weakly by the simultaneous-approximation-term (SAT) 
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technique mmmm- This is often referred to as the SBP-SAT methodology, 
and a summary of it is made in [J. The energy method is used to determine 
the strength of the boundary enforcement and to derive an energy estimate for 
stability. Another candidate for the approximation of spatial derivatives is dis¬ 
continuous Galerkin method, and it has been successfully used in [5] for the 
wave equation. 

In the SBP-SAT framework, accuracy analysis has drawn less attention than 
stability. A commonly used measure for accuracy is convergence rate, typically 
in L 2 norm. The convergence rate indicates how fast the error in the numerical 
solution approaches zero as the mesh size h goes to zero. The error in the 
numerical solution is caused by truncation error. The truncation error near a 
boundary is often larger than that in the interior of the computational domain, 
but with the large boundary truncation error localized at a fixed number of grid 
points that is independent of mesh refinement. As a consequence, its effect to the 
convergence rate may be weakened. It is well-known that by applying directly 
the energy method to the error equation, 1/2 order is gained in convergence. 
However, this convergence result is in many cases suboptinral. When using a 
finite difference method with mesh size h to solve a first order hyperbolic partial 
differential equation, it is in many cases possible to show that (p + l) th order 
convergence rate can be obtained if the boundary truncation error is 0(h p ) 
mm- In other words, we gain one order in convergence rate. The technique 
used to analyze the effect of the boundary truncation error in [B] is normal mode 
analysis. The idea is based on Laplace transforming the error equation in time. 
The error equation due to boundary truncation error can be formulated as a 
system of linear equations, which is referred to as the boundary system. A gain 
of one order in convergence rate follows if the boundary system is nonsingular 
for all Re(s) > 0, where s is the dual variable of time in Laplace space. For such 
cases we use the same terminology as in [8] and say that the boundary system 
satisfies the determinant condition. Also note that if a problem satisfies an 
energy estimate, the boundary system is nonsingular for Re(s) > 0. However, 
an energy estimate does in general not imply that the determinant condition is 
satisfied. 

Also for many other systems the boundary truncation errors have less severe 
effect than predicted by energy estimates. In Q] a parabolic problem is consid¬ 
ered. Numerical experiments show a convergence rate min(2p,p-|- 2), where p 
is the order of the boundary truncation error, while a careful deviation of the 
energy estimate yields min(2p,p-|- 3/2) as the convergence rate. This reaffirms 
that the energy estimate is an upper bound of the error, and indicates that the 
energy estimate is not always sharp. 

A more general result on convergence rate is presented in El, where it is 
shown that k orders can be recovered from boundary truncation error for a 
partial differential equation with k th spatial derivatives. Parabolic, incomplete 
parabolic and second order hyperbolic partial differential equations are consid¬ 
ered, and for these equations we call min(2 p,p + 2) the optimal convergence 
rate. The condition for the optimal convergence is formulated in terms of a new 
stability concept, pointwise stability, but the underlying analysis of the optimal 
convergence is based on the determinant condition, and therefore the approach 
cannot be applied when the boundary system is singular at the origin. 

When using SBP-SAT finite difference method to solve second order hy¬ 
perbolic partial differential equations, there are many cases that the boundary 
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system does not satisfy the determinant condition. In this paper, we aim at fill¬ 
ing the gap in the study of convergence rate for such problems. We use the wave 
equation in second order form as our model problem. In particular, we consider 
problems with Dirichlet boundary condition, Neumann boundary condition and 
grid interface condition. We refer to them as the Dirichlet problem, the Neu¬ 
mann problem and the interface problem, respectively. By analyzing the error 
equation, we obtain a boundary system of linear equations. If the determinant 
condition is satisfied, we would get the optimal convergence rate as shown in 
[5J Ch.12]. If the determinant condition is not satisfied on the imaginary axis, 
mm are not applicable and one has to carefully derive an estimate for the 
solution of the boundary system to see how much is gained in convergence rate. 

In this paper, we find that 1) an energy estimate for the second order wave 
equation does not imply an optimal convergence rate; 2) the determinant con¬ 
dition is not necessary for an optimal convergence rate; 3) if there is an energy 
estimate but the determinant condition is not satisfied, there can be an optimal 
gain of order 2, or a non-optimal gain of order 1, or only a 1/2 order gain that 
is given by the energy estimate. We demonstrate how to analyze the boundary 
system to determine the accuracy gain. 

More specifically, for the Dirichlet problem, the boundary system is singular 
or nonsingular depending on the value of the penalty parameter in SAT. In this 
case, the optimal convergence rate is only obtained when the determinant condi¬ 
tion is satisfied. For the Neumann and the interface problems, the determinant 
condition is not satisfied even though the schemes are stable. We will analyze 
the boundary system, and show how much is gained in convergence in those 
cases. The analysis agrees well with numerical experiments. In [151 [151 , the 
same technique is also used to prove optimal convergence rate for Schrodinger 
equation. 

The structure of this paper is as follows. We start in [J5]with the SBP-SAT 
method applied to the one dimensional wave equation with Dirichlet boundary 
condition. We apply the energy method and normal mode analysis, and derive 
results that correspond to the second, fourth and sixth order schemes. We 
then discuss the Neumann problem in 21 In 21 we consider one dimensional 
wave equation on a grid with a grid interface. We extend the analysis to two 
dimensions in 2) In 21 we perform numerical experiments. The computational 
convergence results support the accuracy analysis. The conclusions are drawn 
in 21 


2 One dimensional wave equation with Dirichlet 
boundary condition 

To describe the properties of the equation and the numerical methods, we 
need the following definitions. Let wi(x ) and W 2 (x) be real-valued functions 
in L 2 [ 01 , 02 ]- The inner product is defined as (uq , 1 x 2 ) = w\w^dx. The 
corresponding norm is ||uq|| 2 = (uq,wq). 

The second order wave equation in one dimension takes the form 


U t t = U xx + F, x G [0, 00 ), t > 0, 
U{x,0) = f(x), U t (x, 0) = f 2 (x). 


( 1 ) 
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where F is a forcing term. We consider the half line problem with Dirichlet 
boundary condition: 

U(0,t) = g(t). (2) 

The forcing function F, the initial data f 1 , f 2 and the boundary data g are 
compatible smooth functions with compact support. For the problem © in a 
bounded domain, there is one boundary condition at each boundary. For the 
half line problem in consideration, the right boundary condition is substituted 
by requiring that the L 2 norm of the solution is bounded, i.e. ||t/(-,f)|| < 00 . 
The problem (fT j) - © is well-posed if there is an energy estimate with g(t) = 0 
and the problem is boundary stable. In in Definition 2.3], it is defined that 
the problem ©-© is boundary stable if with f 1 = f 2 = 0 and F = 0 there are 
constants 770 > 0, I\ > 0 and a > 0 independent of g such that for all 77 > 770 , 
T > 0, 

[ e - 2 *\u(0,t)\ 2 dt<^ f e~ 2pt \g(t)\ 2 dt. 

Jo V Jo 

It is obvious that © is boundary stable with Dirichlet boundary condition. To 
derive an energy estimate, we multiply Equation © by Ut and integrate by 
parts. With homogeneous Dirichlet boundary condition g(t) = 0, we obtain 

fv^ <11711, 

where E = \\Ut\\ 2 + ||bfr|| 2 is the continuous energy. The energy estimate follows 
from Gronwall’s lemma 

V^< v / ll/ill 2 + ll/T+ /Vo,411*. (3) 

JO 

Therefore, problem © © is well -posed. 

Next, we introduce the equidistant grid Xi = ih, i = 0,1,2, , and a 

grid function Ui(t) ~ U(xi,t). Furthermore, let u(t) = [uo(i), ui(t), ■ ■ ■ ] T . We 
also define an inner product and norm for the grid functions a and b in 1Z 
as (a, b)n = a T Hb and \\a\\ 2 H = a T Ha , respectively, where a T denotes the 
transpose of a and H is a positive definite operator in the space of grid functions. 
An SBP operator approximates a derivative, and mimics the property of the 
continuous integration-by-parts via the inner product and norm defined above. 

In this paper, we use the diagonal norm SBP operator D approximating 
second order derivative, i.e. D rj The order of accuracy is denoted by 

2 p, p = 1,2,3,4, 5. For a 2 p th order diagonal norm SBP operator mm, the 
error for the approximation of the spatial derivative is 0(h 2p ) in the interior, and 
0(h p ) near the boundary. The SBP operator can be written as D = H~ 1 (—A + 
BS ), where H is diagonal and positive definite, A is symmetric positive semi- 
definite and B = diag(— 1,0, 0, • • •). S' is a one sided approximation of the first 
derivative at the boundary with the order of accuracy p+ 1. The elements in D 
are proportional to 1 /h 2 while the elements in H~ l and S are proportional to 
1 jh. Another version of the SBP operators are constructed with block norms, 
see |TT. Q . 

We use the SAT technique to impose the Dirichlet boundary condition 
weakly. The semi-discretized equation corresponding to © and © with ho¬ 
mogeneous boundary data reads: 

Utt = Du — U -1 S T Equ — —H~ 1 Equ + F g , 


( 4 ) 
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2 p 

2 

4 

6 

8 

10 

a 2p 

0.4 

0.2508560249 

0.1878715026 

0.0015782259 

0.0351202265 


Table 1: a 2p values 


where eo = [1, 0,0, • • • ] T , Eo = eoe^ and F g is the grid function corresponding to 
F(x, t). On the right hand side, the first term is an approximation of U xx , while 
the second and third terms impose weakly the boundary condition 17(0, t) = 0. 
They act as penalty terms dragging the numerical solution at the boundary 
towards zero so that the boundary condition is simultaneously approximated. In 
general, the boundary condition is not satisfied exactly. The penalty parameter 
t is to be determined so that an energy estimate of the discrete solution exists, 
which ensures stability. 


2.1 Stability 

In [2 [121 13], it is shown that the operator A can be written as 

A = ha 2p (E 0 S) T E 0 S + A, (5) 


where a 2p > 0 is as large as possible and A is symmetric positive semi-definite. 
This is often called the borrowing trick. We use the technique in m to compute 
the values of a 2p and list the results in Table [T] 

Multiplying Equation (JJ) by uf H from the left, with (JSJ) we obtain 


/ 


d_ 

dt 


\ 


\ u t\\ 2 H + W U W\ + yV^ a 2 P (BSu) 0 — 


U 0 


Ot-2p 


V 


= 2 u t HF a 


E h 




For t > A-, we have Eh > 0. In this case, |||u|]|^ := Eh is a discrete energy, and 
Hl'lll^ is the corresponding discrete energy norm. We then obtain the discrete 
energy estimate 

< \JE h ,o + [ \\F g (-,z)\\ H dz, 

Jo 

where Eh,o is the initial discrete energy. Clearly, stability is ensured if r > 

r 2p a 2 p ' 

Remark It is possible to include the term ||ii||^ in the discrete energy Eh- A 
similar energy estimate holds [8]. This is useful when comparing the numerical 
results with the theoretical analysis. 



2.2 Accuracy analysis by the energy estimate 

Assuming that the continuous problem has a smooth solution U(x,t), we can 
derive the error equation for the pointwise error = U(xj , t) — Uj (t ). The error 
equation corresponding to (0} is 

& = D£- H~ 1 S t E 0 ^ - Ih-'EoZ + T 2p , 


( 6 ) 
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where ||£(f)||/ l < oo, T 2p is the truncation error and 2 p is the accuracy order of 
the SBP operator. The first m components of T 2p are of order 0(h p ), and all the 
other components are of order 0(h 2p ). For 2 p = 2,4, 6,8,10, the corresponding 
m values are 1,4,6, 8,11. In the analysis, we only consider the leading term of 
the truncation error. We introduce the interior truncation error T 2p,/ , and the 
boundary truncation error T 2p ’ B such that T 2p = h 2p T 2p ’ 1 + h p T 2p ’ B , where 
T 2p ’ T and T 2p ' B are independent of h, but depend on the derivatives of U. We 
have 

t 2 P ,i\ = 0 0<*<m-l t 2 P ,b |=0 i>m- 1 

1 0 i>m— 1 1 0 0<*<m—1 

Since the number of grid points with the large truncation error 0(h p ) is finite 
and independent of h, we have 

oo m— 1 

||T 2p ||/ l = h{h 4p Y \ T i P,I \ 2 + h2p Y I T i P ’ B \ 2 ) < I<ih 4p + K B h 2p+1 . 

i=m i=0 

Si v ✓ S v ✓ 

Oih- 1 ) 0(1) 

For 2 p = 2,4,6, 8,10 and small h, the first term is much smaller than the second 
one. Thus, we have ||T 2p ||? l < K B h p+1 / 2 . By applying energy method to the 
error equation ©, we obtain 

m\ h < c e h p ^' 2 . 

This means that by the energy method we can only prove a gain in accuracy 
order of 1/2. Therefore, the convergence rate is at least p + i if the numerical 
scheme is stable, that is if r > t- 2 P . 

2.3 Normal mode analysis for the boundary truncation 
error 

To derive a sharp estimate, we partition the pointwise error into two parts, the 
interior error e 1 and the boundary error e, such that £ = e 1 + e. e 1 is the error 
due to the interior truncation error, and e is the error due to the boundary 
truncation error, e 1 can be estimated by the energy method, yielding 

infill h <C!h 2p . (7) 

where C/ is a constant independent of h. 

We perform a normal mode analysis to analyze the boundary truncation 
error e for second, fourth and sixth order SBP-SAT scheme. The boundary 
error equation is 


e a = De- H~ 1 S T E 0 e - 1-H~ l E 0 e + h p T 2p ’ B , (8) 

h 

where ||e||h < oo. In the analysis of e, we take the Laplace transform in time of 

m, 

s 2 e = De - H~ 1 S T E 0 e - yH^Eoe + h p T 2p ’ B , 

a 


(9) 
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2 p 

q(,r > r 2p ) 

q(T = t 2 p) 

2 

2 

1.5 

4 

4 

2.5 

6 

5 

3.5 


Table 2: Theoretical convergence result for one dimensional wave equation. 


for Re(s) > 0, where' denotes the variable in the Laplace space. Let s = sh. 
After multipling h 2 on both sides of ©, we obtain 

s 2 i = h 2 De - h 2 H~ 1 S T E 0 e - ThH~ x E 0 e + h p+2 T 2p ’ B . (10) 

Note that the first three terms in the right hand side of m are /i-independent. 

By assumption the true solution U (x, t ) is smooth and, consequently, | T 2p,B ( s ) | 
decreases fast for large |s|. Convergence rate is an asymptotic behaviour as the 
mesh size h approaches zero. Therefore, we need to investigate the solution of 
(USD as s approaches zero. In the analysis, we assume |s| < K and consider 
Re(s) > ?? > 0, where K and are some positive constants. Equivalently, 
|s| < 1 and Re(s) > rjh. If the semi-discretized equation is stable, the Laplace 
transformed problem is nonsingular for Re(s) > 0 [5]. 

There are essentially two steps in the normal mode analysis. Firstly, by a 
detailed analysis of the error equation m, we obtain an estimate for ||e||h in 
the Laplace space. Then we use Parseval’s relation to derive an estimate for the 
error in the physical space of the form 



\\e(;t)\\ 2 h dt<Ch\ 


( 11 ) 


where C depends only on 77 , the final time tf and the derivatives of the true 
solution U. The results for the Dirichlet problem are summarized in Theorem 

nm 


Theorem 2.1. For the second, fourth and sixth order stable SBP-SAT approx¬ 
imations of the second order wave equation am on a half line with Dirichlet 
boundary condition, the rates q in ill)) depend on t, and are listed in Tabled 

In )}51 we present results from numerical experiments, which agree well with 
the theoretical results in Table [2] After some preliminaries we will prove The¬ 
orem 12.11 for 2 p = 2,4 and 6 in T2.3.1I 12.3.21 and 12.3.31 respectively. In the 
proof, we will show that the determinant condition is satisfied if r > r^p but 
not satisfied if r = T 2 p. In addition, the determinant condition is necessary for 
the optimal convergence rate. 

To begin with, we note that sufficiently far away from the boundary, the 
two penalty terms and the boundary truncation error in m are not present. 
The coefficients in D correspond to that in the standard central finite difference 












scheme. We have 


2p = 2 : s 2 ej = D + D_£j , j = 3,4, 5, • • •, (12a) 

2p = 4: s 2 e j = (D + D_--^(D + D_) 2 )e j , j = 4,5,6,---, (12b) 

2p = 6: s 2 ij = (D + D_ - ^(. D+D.) 2 + ^( D+D _) 3 )e J -, j = 5, 6, 7, • • •, 

(12c) 

■where 

an d D _e j = kzlt±. 
h h 

The corresponding characteristic equations are 

2p = 2 : k 2 — (2 + s 2 )k + 1 = 0, (13a) 

2p = 4 : k 4 — 16k 3 + (30 + 12S 2 )k 2 - 16k + 1 = 0, (13b) 

2p = 6 : 2k 6 - 27k 5 + 270k 4 - (180s 2 + 490 )k 3 + 270k 2 - 27k + 2 = 0. 

(13c) 

It is easy to verify by von Neumann analysis that the interior numerical scheme 
is stable when applied to the corresponding periodic problem. From Lemma 
12.1.3 in (8) pp. 379], it is straightforward to prove that there is no root with 
|k| = 1 for Re(s) > 0. We will need the following specifics for the roots. 

Lemma 2.2. For 2p = 2,4,6, the number of admissible roots of m satisfying 
|k| < 1 for Re(s) > 0 is 1,2,3, respectively. In the vicinity of s = 0, they are 
given by 

2p = 2: ki = l-s + C>(s 2 ), (14a) 

2p = 4: ki = 1 - s + 0(s 2 ), k 2 = 7-4V3 + C>(s 2 ), (14b) 

2p = 6: ki = 1 - s + 0(s 2 ), k 2 = 0.0519- 0.0801* + £>(s 2 ), (14c) 

k 3 = 0.0519 + 0.0801* + 0( s 2 ). 

Proof. 2p = 2: Equation (113al) has two roots: ki ,2 = 1 + ^s 2 ± ^Vs 4 + 4s 2 . We 
find by Taylor expansion at s = 0 that Ki = 1 — s + 0(s 2 ) and k 2 = l + s + C>(s 2 ). 
Thus, the admissible root is n\ = 1 — s + 0(3 2 ). 

2p = 4: Equation (I13bl> has four roots: 

ki = \J2A- 3s — 8\/9- 3s 2 - \/9 - 3s 2 + 4, 
k 2 = 4 — a/ 9 - 3s 2 - ^24- 3s- 8\/9- 3s 2 , 
k 3 = V9 - 3s 2 - \J 8\/9 — 3s 2 — 3s + 24 + 4, 
k 4 = \J 8\/ 9 — 3s 2 — 3s + 24 + y/9 - 3s 2 + 4. 

We find by Taylor expansion at s = 0 that 

Ki = 1 - s + 0(s 2 ), k 2 = 7-4v / 3 + £>(s 2 ), 

k 3 = 7 + 4 v / 3 + 0(s 2 ), k 4 = 1 + s + 0(s 2 ). 
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Thus, the admissible roots are m = 1 — s + 0(s 2 ) and K 2 = 7 — 4y/3 + 0(s 2 ). 

2p = 6: (I13cl) is a six order equation. There is no formula for general six 
order equations over the rationals in terms of radicals. We solve numerically 
equation itrsa with s = 0, and then analyze the six roots by perturbation 
theory to find out their dependence on s. The three admissible roots are given 
by (114cl) . □ 

Lemma 2.3. Consider |s| -C 1 and Re(s ) > ph > 0. For 2 p = 2,4,6 the 
admissible roots Ki(s) satisfy 

1 — |«i(s)| 2 > 2Re[s) + 0(|s| 2 ). 

Proof. Let s = x + iy where x, y are real numbers. Then \x\ > ph. 
l-| Kl (s)| 2 = l-|l-s + 0(S 2 )| 2 

> 1 _| 1_ 5 |2 + 0 ( 15 ( 2 ) 

= l-\l-x-yi\ 2 + 0(\s\ 2 ) 

= 2 x — x 2 — y 2 + 0(|s| 2 ) 

= 2Re(s) + 0(|s| 2 ). 


□ 

By Lemma ROl we obtain 1 _| k 1 ^^ 2 < ^ to the leading order. We will use 
this inequality to estimate ||e||/j. 

2.3.1 Proof of Theorem I2TT1 for the second order scheme 

Away from the boundary, equation (flOll is simplified to (112 a I) . which is a recur¬ 
sion relation for ij. A general solution of USD can be written as 

ej =an{~ 2 , j = 2,3,4,-, (15) 

where Ki is obtained from (I14al) and a is an unknown in the boundary system. 
In the following, we derive the boundary system and investigate its solution 
near s = 0 to obtain an estimate of the error in the Laplace space. The first 
three rows of ohd are affected by the penalty terms. They are 

s 2 eo = £o — 2ei + £2 + 3fo — 2r£o + h 3 T 0 ' , 

S 2 £i = £o — 2 £i + £2 — 2fio, 

S 2 £2 = £l - 2 e 2 + £3 — -£o- 

By Taylor expansion, it is straightforward to derive Tq' B = — U X xx(S>t s) to the 
leading order. We write m in the matrix vector multiplication form with the 
help of (fl5ll . and obtain the boundary system 


— 1 s 2 — 4 + 2r 

-1 1 

Ki — 2 — s 2 \ 

C/31 

to 

I -1 _|_ tfO 

to 

1_1 


a 

£o 

_£i_ 

V 

II 

w 

f 

0 

0 

C 2 d(s,t) 



^2 £) 


r£i 2 , B 
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Next, we investigate the invertibility of C 2 d{s , r) in a vicinity of s = 0 to derive 
an estimate for |S 2 _d|, where | • | denotse the standard Euclidean norm of vectors 
and matrices. We note that 


C 2 d( 0 ,t) 


-1 

-1 

-1 


-4+2 r 2 
1 2 


The determinant of C 2 ,d( 0 ,t) is given by det(C 2 ,o(0, r)) = 5 — 2r. Clearly, 
is singular if r = 2.5 and nonsingular otherwise. The stability condi¬ 
tion is r > t 2 = 2.5. 

In the case r > r 2) C 2 d( 0 ,t) is nonsingular. By the continuity of C 2 d(s, t), 
for any 7 > 1 and sufficiently small h we have |C^(s, r)| < 7 |C^( 0 , r)|. Thus, 
|S 2 d| < 7|C“^(0,t)||T^’ B |/i 3 . Therefore, we have \a\ 2 < K a \U xxx (0, s)\ 2 h 6 and 
E- = o N 2 — Ke\U xxx (0, s)| 2 h 6 , where K a and K e are constants independent of 
h. In L 2 norm, we have 


llell 2 = h^2\ei\ 2 + h\a\ 2 ]T \m\ 21 = h N 2 + h W\ 


2=0 


2=0 


2=0 


1 - \Ki 


12 • 


II B,h 


By Lemma [231 we have ||e||f j?l = h\a \ 2 < ^\U XXX (0, s)\ 2 h 6 . The first 

term can be bounded as ||e||| h = hSi= 0 N 2 < K e \U xxx (0, s)\ 2 h 7 < ||e||f ^ for 
small h. Thus, to the leading order, 

\\e\\ 2 h <^\U xxx (0,s)\ 2 h e , 

where K a is a constant independent of h. By Parseval’s relation, we have 

/»oo -1 /»oo 

J o e ~ 2rlt Mhdt = — J P(?7 + *0II hdZ 

\U xxx (0,r] + i£)\ 2 di 


< 


' —OO 

k„h 6 r°° 

J-00 

k a h e r°° 
Jo 


e~ 2r,t \U xxx (0, t)\ 2 dt. 


Note that in the above, we use Parseval’s relation twice. By arguing that the 
future cannot affect the past, we can replace the upper limit of the integrals on 
both sides by a finite time tf. Since 

[ tf e~^ tf \\e\\ 2 dt < /V^Uell 2 dt, 

Jo Jo 

[ ' e~ 2rit \U xxx (0,t)\ 2 dt < [ ' \U xxx (0,t)\ 2 dt, 


we obtain 


\e\\ 2 h dt < h 3 \ 


I k a e 2vt f r tf 

Jo 


\u xxx {0,t)\ 2 dt. 
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Thus, the boundary error is 0(h 3 ). In this case, the interior error is 0(h 2 ) by 
0. which is the dominating source of error. Therefore, the overall convergence 
rate is 2 . 

In the case r = t 2 , C 2 B (0, t) is singular, and the above analysis is not valid. 
By the energy estimate the convergence rate is p + ^ = 1.5. 


(17) 


2.3.2 Proof of Theorem I2TD for the fourth order scheme 

Away from the boundary, the solution of (fTTfll is 

e? = crin [- 2 + a 2 K J 2 ~ 2 , j = 2,3,4, • • • , 

where n \, k 2 are obtained from (I14bll , and er 4 , ct 2 will be determined by the 
boundary system. The first four rows of dll are affected by the penalty terms. 
They are 

/122 — 48t' 


s 2 e 0 = 


V 17 

85 


s ei — ——£q — 2 ei + e 2 + h T j 


Co — 5ei + 4e 2 — 63 + /i 4 Tq’ B , 




59 

_ 2 . 68 . 59. 

S £2 = 43 C ° + 43* 


17, 


59. 


Se3 = ~49 eo+ 49 e2 


110. 59. 4 4 . 4)B 

-69 -b —63 — —64 H - ri 1 9 , 

43 43 43 25 

118. 64. 4 4 - 4)B 

dT 3 + 49* “ 49 65 + h T * ■ 


(18) 


By Taylor expansion, we have 

n4 ,B 


Tq’" = hj7 rara ( 0,s), T 4 ’ s = --Lt/ xxra (0, s ), 

n4,S 


h4,B 


L 2 — 516 ^ xxxx ^’’ 5 )’ ^3 — ^gg U XX xx (O 7 5 ) 

to the leading order. With s = 0, the boundary system follows from 03 OHD 

C 4B (0,t)£ 4B = h 4 T 4 ’ B , 

where £ 4B = [< 7 i, ct 2 , e 0 , e 4 ] T , T 4,B = [Tq ,b ,T 4 ’ B ,T 2 ’ B ,f^’ B ] T and C 4 d( 0 ,t) 
is given in Appendix 1 . We solve for det(C 4 B (0, r)) = 0, and obtain r = 
2 77 ^g^^T) 9 ^ ~ 3.986350342. Recall that the stability condition is r > r 4 ~ 
3.986350342. 

In the case r > r 4 , C 4 B (0,t) is nonsingular. Similarly to the second order 
case in 1 )2.3.11 we obtain to leading order 

||e||^<^|^ ra ,( 0 , S )| 2 / l 8 . 

By Parseval’s relation and by arguing that the future cannot affect the past, we 
obtain 


J \\e\\ 2 h dt<h 4 ^ 


I K ai e 2??t / r*t 

Jo 


\U xxxx (0,t)\ 2 dt. 


Thus, the boundary error is 0(h 4 ). In this case, the interior error is also 0[h 4 ) 
given by ©. Therefore, the convergence rate is 4. 

In the case r = r 4 , C^i^O, r 4 ) is singular. By the energy estimate, the 
convergence rate is p + 1 


1 =2.5. 
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2.3.3 Proof of Theorem 12.11 for the sixth order scheme 

Away from the boundary, the solution of (1101) is 

e? = CTi/tf 3 + ct 2 a 4 -3 + ct 3 k| -3 , j = 3,4,5, •• • , (19) 

where Ki,k 2 ,k 3 are obtained from (I14cl) . and <7 i,<j 2 ,<7 3 are determined by the 
boundary system. Near the boundary, the first six rows of m are affected by 
the penalty terms. In the same way as for the second and fourth order method, 
we obtain the boundary system in = [eo, ei, e 2 , <7i, <r 2 , cr 3 ] T : 

C 6D (s,T)X 6D = h 5 T^ B , 

where Cqd{s,t) is given in Appendix 1. We find that C6 d(0,t) is singular if 
r = 76, and nonsingular otherwise. When r > tq, we use Taylor expansion to 
express T® ,B in terms of the derivatives of the true solution £/, and Lemma I2T51 
to derive an estimate for (HU). The boundary error is 0(h 5 ). In this case, the 
interior error is 0(h e ) given by ([7]). Thus, the convergence rate is 5. 

When t = tq, the convergence rate p + ^ = 3.5 is given by the energy 
estimate. 

3 One dimensional wave equation with Neumann 
boundary condition 

Next, we consider an example which never satisfies the determinant condition, 
but nonetheless exhibits optimal convergence. We consider Equation m with 
Neumann boundary condition 


£4(0, t) = g(t). (20) 

We use the same assumption of the data as for the Dirichlet problem. With 
homogeneous Neumann boundary condition g(t) = 0 we get an energy estimate 
identical to the one for the Dirichlet problem |3]). To prove well-posedness, we 
also need to show that the equation is boundary stable. Boundary stability 
for the Neumann problem is proved in m Theorem 3.8] by Laplace transform 
technique. 

3.1 Stability 

To discretize the equation in space, we use the grid and grid functions intro¬ 
duced for the Dirichlet problem. The semi-discretized equation of the Neumann 
problem is 

u tt = Du + H^EoSu + Fg. (21) 

On the right hand side of CD). the first term approximates U xx and the second 
term imposes weakly the boundary condition U x (0,t) = 0. Since we consider a 
half line problem, the SBP operator D can be written as D = H~ 1 (—M — EoS). 
As a consequence, the semi-discretization CD can be written as 

u tt = -H~ 1 Mu + F g . 


( 22 ) 
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Multiplying Equation (12^1) by uf H from the left, we obtain 
j t {\\u t \\ 2 H + \\u\\ 2 M ) = 2u t HF. 

The discrete energy E^ eu = |||'tt|||^ eli := ||wt|||f + \\u\\m * s bounded as 



m is the energy estimate. 

3.2 Accuracy 

Define the pointwise error as £, = U(xj,t) — Uj(t). The error equation corre¬ 
sponding to (l22l) is 

£ tt = -H~ 1 M£ + T 2p ’ Neu , (24) 

where ||£||h, < oo, T 2p,Neu is the truncation error and 2 p is the accuracy order 
of the SBP operator. In the same way as for the Dirichlet problem, by applying 
the energy method to the error equation (EH) , we obtain an estimate 

Klllf “ < C Neu h p+1 

where C^eu is a constant independent of h. This means that the convergence 
rate of (El is at least p + 1/2. In the following, we use normal mode analysis 
to derive a sharp bound of the error, which agrees well with the results from 
numerical experiments. 

Similar to the Dirichlet problem, we partition the error £ into two parts, 
the error e 1 due to the interior truncation error and the error e due to the 
boundary truncation error, e 1 can be estimated by the energy method, yielding 

lll^mNeu < fjNeu^ 2 p w j iere (jNeu a cons t arL t independent of h. The boundary 

error equation is 

e tt = -H~ 1 Me+h p T 2p ’ Neu ’ B (25) 

where hPT 2p ’ Neu ’ B is the boundary truncation error. T 2p,Neu ’ B depends on 
the derivatives of the exact solution U, but not h. Moreover, only the first m 
elements of T 2p ’ Neu ’ B are nonzero, where m depends on p but not h. 

We Laplace transform Equation (l25l) . with the notation s = sh , we obtain 

s 2 e = —h 2 H~ 1 Mi + h p+2 T 2p ’ Neu ’ B (26) 

Comparing with the error equation for the Dirichlet problem (HOD , there are 
less rows affected by the boundary closure than that for the Drichlet problem. 
The characteristic equations of the interior error equations are the same as 
for the Dirichlet problem m- Lemma El and EH are also applicable to the 
Neumann problem. For the error equation near the boundary, we analyze below 
the second, fourth and sixth order accurate cases. 




14 


3.2.1 Second order accurate scheme 

For second order accurate scheme, only the first row of the error equation is 
affected by the boundary closure. The general solution of the error equation is 

e j = ( TK{,j= 0,1,2, - - - , (27) 

where = 1 — s + 0(s 2 ) is obtained from (114aD and a will be determined by 
the boundary system. The first row of (EHD is 

s 6o = —2e 0 + 2ei + h°T 0 

From Equation (ED), we obtain 

C 2N {s)a = h 3 T 2 ’ Neu ' B , C 2N (s ) = s 2 + 2 - 2 Kl = 23 + 0(s 2 ). (28) 

Clearly, C 2 jv(0) = 0. By Taylor expansion, we have 

T o’ NeU ’ B = ~U xxx (0,s) + -U xxx (0,s) = --U xxx (0,s). 

where —h 3 U xxx { 0, s ) is the truncation error by the SBP operator and h 3 U xxx (0 , s)/3 
is the truncation error by the SAT. Therefore, we have |<r| < \U XXX (0, s)\h 2 . 

In the same way as for the Dirichlet problem, an estimate in physical space 
follows, and we find that the boundary error is 0(h 2 ). Since the interior error 
is also 0(h 2 ), the convergence rate is 2. For second order accurate scheme, we 
gain two orders from the boundary truncation error in the Dirichlet problem, 
but only gain one order in the Neumann problem. Fortunately, this does not 
affect the overall convergence rate. 

3.2.2 Fourth and sixth order accurate scheme 

We follow the above approach to derive the four by four boundary system 

CW(S)E 4 jv = h*fZ’ Neu ' B , 

for the fourth order scheme. At s = 0 we find that C4vr(0) is singular with rank 
3. By Taylor expansion, we derive the precise form of ’jr^’ Neu,B . To derive a 
sharp estimate for |£ 4 at|, we write C 4 jv(s) = C 4 at( 0 )+^( 74 ^( 5 ) + C?(s 2 ) and use 
Lemma 3.4 in [13] . Let singular value decomposition of C4jv(0) be UI n S±nV±n■ 

We find that Ul N C' iN (Q)ViN = —0.3095 ^ 0. Moreover, T^’ Neu ’ B is in the 
column space of C4jv(0). Hence, |E 4 at| < K4n\U xxxx ( 0, s)\h 4 where AAat is a 
constant independent of h. This means that the boundary error is (7(/i 4 ) and 
there is a gain of 2 in accuracy order. In this case, the interior error is also 
0(h 4 ), and the overall convergence rate is 4. We show C4jv(0), 6^(0) and 
f4,Neu,B in Appendix 2. 

The situation of the sixth order scheme is very similar to that of the fourth 
order scheme. The six by six boundary system 

CW(S)E ( sjv = h 5 f%’ Neu ’ B 

is singular at s = 0 with rank 5. The exact form of T®' NeU)B is derived by 
Taylor expansion. Here again we find that Lemma 3.4 in H5I is applicable with 
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^6N^6n(.^)^n = 0-2072 — O.OOOli 7 ^ 0, and f^’ Neu ’ B is in the column space of 
C 6 iv( 0 ). As a consequence, the boundary error is 0(h 5 ). Since the interior error 
is 0(h 6 ), the convergence rate is 5. C6jv(0), C' 6N (0) and T® ,JVeu,B are given in 
Appendix 2. 

Remark For fourth and sixth order scheme, the optimal p + 2 convergence rate 
relies on the fact that the boundary truncation error vector is in the column 
space of C 4 / 6 iv( 0 ). In the numerical experiments in we will verify the optimal 
convergence rate. If the boundary truncation error would not be in the column 
space of C 4 / 6 jv( 0 ), we only obtain p + 1 convergence rate by using Lemma 3.4 
in m■ In the numerical experiments, we will add a dissipative term to the 
boundary part of the SBP operator, so that the boundary truncation error is 
not in the column space of CW^O) but remains the same magnitude. We 
indeed observe (p + l) th order convergence rate. 


4 One dimensional wave equation with a grid 
interface 


In this section, we consider another example that does not satisfy the determi¬ 
nant condition. It is the Cauchy problem for the second order wave equation in 
one dimension 


U t t = U xx + F, x G (- 00 , 00 ), t> 0, ||J/(-,f)|| < 00 
U(x,0) = f\x), U t (x, 0) = f 2 {x), 


where the forcing function F , the initial data / 1 and / 2 are compatible smooth 
functions with compact support. It is straightforward to derive the energy 
estimate of the form ©• 

We solve the equation on a grid with a grid interface at x = 0. With 
the assumption that the exact solution is smooth, it is natural to impose two 
interface conditions at x = 0: continuity of the solution and continuity of the 
first derivative of the solution. 

We introduce the grid on the left X-j = — jhi, , j = 0,1, 2, • • •, and the grid 
on the right Xj = jh,R, j = 0,1, 2, • • •, with the grid size h R and Hr, respectively. 
The grid interface is located at x = 0. The grid functions are u B (t) ~ U(x-j,t ) 
and uf{t) ~ U(xj,t). Denote 


eo l = [■■•, 0, 0,1] T , e 0 R = [1, 0,0, • • • ] T , 

U L = [ ,u£ 1 ,uS'] T , U* = [«?,«?,-f- 


Both Uq = ej l u l and = e B R u R approximate 17(0, t). Both (Slu l ) 0 = 
e oand ( Sru r )q = e R R Snu R approximate U x (0,t). To simplify notation, 
we define {«} = Uq —u r and {Su} = (Slu l ) 0 — ( Sru r )q. The semi-discretized 
equation reads 

Uu = d lu l + e 0 L {u} - ^H^eorASu} - ^o^u} + F gL , 

u tt = Dru r + -H R 1 S r cor{u} — -H^eoRiSu} — TRH^eoRiu} + F gR , 

(30) 
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2 p q(r > f 2p ) q(r = f 2p ) 

2 2 1.5 

4 4 2.5 


Table 3: Theoretical convergence result for one dimensional wave equation with 
a grid interface. 


where FgL and F gR are the grid functions corresponding to F(x , t). The penalty 
parameters tl and tr are chosen so that the semi-discretization is stable. By the 
energy method, the numerical scheme is stable if tl = —tr > ^ZtVZjiR '■= f 2p . 
The stability proof is found in [12] Lemma 2.4, pp. 215]. If Jir = Hr := h , then 
T2 P = 2a 2p h- 

By the energy estimate, the convergence rate is at least p + \ if the semi¬ 
discretization is stable. In order to derive a sharper estimate, we follow the 
same approach as for the half line problem in In the remaining part of this 
section we will only consider the effect of the interface truncation error. The 
interior truncation error results in an error 0(h 2p ) in the solution. Denote 

= U(x-j , t) - uZjft), ef(t) = U(xj , t) - j = 0,1,2,- , 

and t = tr = —tr. The error equation reads 

e tt = D L e L + ^H^SZeoLie} - ^H^eo^Se} - tH~ l \ 0L {e} + h p L T 2p ’ L ’ B , 

4 = d R£ R + \HR l Sle 0R {e} - l -H~ l e 0R {Se} + rH~ l e 0 fl {e} + h p R T 2p ’ R ’ B , 

with the truncation error h p L T 2p ’ L ’ B and h p R T 2p ’ R F as the forcing terms. We 
take the Laplace transform in time of the interface error equation, and obtain 

~s 2 L e L = h 2 L D L e L + ^H-'SZeorAe} - ^H^e^Se} (31) 

- rhlH^eoLie) + h p L +2 T 2p ’ L ’ B , 

~s 2 r£ R = h 2 R D R e R + ^H R l S T R e 0 R{e} - ^H^eo R {Se} (32) 
+ rh 2 R H R l e 0R {e} + h p + 2 T 2p ’ R ’ B , 

where §l = sh R and sr = sh.R. We use normal mode analysis to derive the error 
estimates for the second and fourth order method. The result is summarized in 
the following theorem. 

Theorem 4.1. For the second and fourth order stable SBP-SAT approximation 
w of the Cauchy problem with a grid interface, the numerical solution con¬ 
verges to the true solution at rate q. The values of q depend on t, and are listed 
in Tabled 

We use the notation h = h^ = rhR in the proof where r is the mesh size 
ratio. The proof follows in the same approach as before. We construct the char¬ 
acteristic equation on each side and solve them to obtain the general solution, 
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which are given by m and cm. respectively. The boundary system is for¬ 
mulated by inserting the general solution to the error equation. The accuracy 
order is determined by how the solution of this system behaves with respect to 
h in the vicinity of s = 0. 


4.1 Proof of Theorem 14.11 for the second order scheme 


Near the interface, the last three rows of m and the first three rows of (1321) 
are affected by the penalty terms. The six by six boundary system takes the 
form 

C 2 i(s,r)^ = h 3 f 2 ’ B , 

where 

E = [a 1L ,a 1R , Cq , e^, eg, ef] T , 

= [§ + 3 ^ 2 , °, °, -|: - l, 0 , 0} T U XXX (0, s). 

C 2 /(0,r) is singular with rank 5 if r > r 2 , and singular with rank 4 if r = r 2 . 
The Taylor expansion of C 2 /(s,r) at s = 0 is given in Appendix 3. 

In the case r > f 2 , we use the Taylor expansion of C 2 /(s,r) at s = 0, and 
Lemma 3.4 in [ 151 to derive an estimate for |E|. We need to perform singular 
value decomposition (SVD) on C 2 /(0,r), which is very difficult due to the vari¬ 
able r in the matrix. Instead, we compute numerically the SVD of C 2 /(0, r) for 
r = l.lf 2 ,1.2f 2 , • • • , 5f 2 . In all cases, we find that (U*C' 2I (0, t)V)§q ^ 0, and 
that T 2,B is not in the column space of C 2 j(0,r). Therefore, in a vicinity of 
|s| = 0, we have \C 2 ^(s,t)\ < Kc\s\~ x . Thus, 

|£| < \C 2I x {s,T)\\fZ’ B \h 3 < K c \s\- x \U xxx (0,s)\h 3 < — \U xxx (0,s)\h 2 , 

V 

where Re(s) > hi] > 0 and K c is independent of h and s. For the error e in L 2 
norm \\e\\ 2 h = ||e L ||^ + \\e R \\ 2 hR , we have 


Pll h — h l ^ l^-il 2 + h-R |ef | 2 + h L \cr L \ 2 ^ \r\l \ 21 + h R \a R \ 2 ^ |«i.r| 


2=0 

1 


2=0 

1 


2=0 2=0 
K c i 


2=0 

R\2 _ 

1-|kil| 2 1 — \rir\ 2 


2=0 


h R \<jR | 


< 

ij 2 

K c i 

< —Y 


|^xxx(0,s)| 2 /l 5 + ^|Pxxx(0,s)| 2 /l 4 


//" 


-\U xxx (0,s)\ 2 h\ 


As before, we obtain 


\JMldt < fe 2 ^/ Acl2 3 - \U xxx (0,t)\ 2 dt. 

Thus, the interface error is 0(h 2 ) and the accuracy gain is only one order. In 
this case, the interior error can be estimated to 0(h 2 ) by an energy estimate, 
and the overall convergence rate is therefore 2. 

In the case r = f 2 , the numerical scheme is still stable but the above analysis 
is no longer valid. By the energy estimate, the convergence rate is p+ - = 1.5. 
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4.2 Proof of Theorem 14.11 for the fourth order scheme 

In this case, the eight by eight boundary system takes the form 

C 4 /(s, t)E = h ir f^ B , 


where 

£ = Wil, cr 2 L,criR, cr 2 R, ef] T , 

115 6 1 5 11 115 6 1 5 11 T - r 

204 ~~ 17r 3 ’ _ 12’ 516’ 588’ 204r 4 ~~ 516r 4 ’ 588r 4 ^ xxxx 

The boundary system is also singular in this case. The Taylor expansion of 
C±i{s, t) at s = 0 is given in Appendix 3. We find that 64 /( 0 , r) is singular 
with rank 7 if r > T 4 , and singular with rank 6 if r = t\. 

In the case r > T 4 , Lemma 3.4 in [15] is applicable since T^ B is in the 
column space of 64 /( 0 , r). We have 

|E| < K s \U xxxx (0,s)\h A , 

where Ky, is independent of h and s. As for the second order case, it then 
follows that the interface error is 0(h 4 ), that is, the gain in accuracy order at 
the interface is 2 . In this case, the interior error is 0(h 4 ) given by an energy 
estimate. Therefore, the convergence rate is 4. 

In the case r = T 4 , the numerical scheme is still stable, and the energy 
estimate gives a sharp estimate. The convergence rate is p+ \ = 2.5. 

5 Two dimensional wave equation 

In this section, we consider the second order wave equation in two dimensions: 

U t t = U xx + U yy + F, x G [0, 00 ), y e (- 00 , 00 ), t > 0, 

U(x,y,0) = f x {x,y), U t {x,y,0) = f 2 (x,y). 

We consider the half plane problem in x with Dirichlet boundary condition 

U(0,y,t) = g(y,t). (34) 

The domain in y is the whole real line. For convenience, we require that all data 
are 27r-periodic in y so that we only need to deal with a finite interval i/ £ [0, 27 t] . 
In addition, we assume that the functions F , f 1 , f 2 and g are compatible smooth 
functions with compact support. The L 2 norm of the solution is bounded, i.e. 

11 U 11 < 00 . Similar to the wave equation in one dimension, with homogeneous 
Dirichlet boundary condition the continuous energy E = \\U t \\ 2 + ||f7 x || 2 + 1111 2 
is bounded by the data. 

Next, we introduce the grid and grid functions. In order to simplify notation, 
we assume that the mesh size is h in both x and y direction. The grid is 

Xi = ih , * = 0,1,2,--- , yj = jh, j = 0,1, 2, • • • , N v , 

such that h = and the grid functions are «y(t) ~ U(xi 7 yj,t). We arrange 
the grid functions Uij(t) column wise in a vector u. The general notation used 




(0,s). 









19 


to expand an operator from one dimension to two dimensions is A x = A x ® I y , 
where ® is Kronecker product and I y is the identity operator. With homoge¬ 
neous boundary condition, the semi-discretized equation reads 

Utt = D XX U + DyypU - H~ 1 S T E 0 U - jH-'Equ + F g . (35) 

o 2 

D xx is the standard SBP operator approximating and D yyp is the standard 

o2 

central finite difference operator approximating . In the same way as for 
the wave equation in one dimension, a discrete energy estimate is obtained if 
the penalty parameter r is chosen so that r > f-, where the values of oli p are 
listed in Table Q] 


5.1 Accuracy analysis 

We analyze the accuracy of the numerical scheme by normal mode analysis. 
Denote e the pointwise error due to the boundary truncation error. Then the 
boundary error equation is 

ett = D xx e + D yvp e - H~ 1 S T E 0 e - ^H- x E 0 e + h p T 2p , 

where ||e||/i < oo. We take the Laplace transform in t and Fourier transform in 
y , and obtain for every uj 

s\l u = h 2 D x J u - h 2 H- 1 S T E 0 ~e u - rhH^E^ + h p+2 f 2p , (36) 


where s+ = \Js 2 — h 2 D“ yp , s = sh and D yyp is the Fourier transform of D yyp . 
For different orders of accuracy, we have 


D 


UJ 

yyp 


sin 2 if 2p = 2, 

< sin 2 if (1 +I sin 2 if) if 2p = 4, 

{ sin 2 if (1 + i sin 2 if + ^ sin 4 f) if 2 p = 6. 


Note that D yyp is nonpositive. 

Lemma 5.1. For Re(s) > 0 and D yyp defined above, it follows that Re(s+) > 0. 

Proof. By the definition of principal square root, Re(s+) > 0. Assume Re(s+) = 
0, then s 2 — h 2 D yyp is real and nonpositive, i.e. Re(s 2 — h 2 D yyp ) < 0. We then 
obtain Re(s 2 ) < R e(h 2 D yyp ) < 0. This contradicts Re(s) > 0. Therefore, 
Re(s+) >0. □ 

The truncation error T 2p can be written as T 2p = [T f 2p ; T 2p ;•••], where 
T 2p = [T 2p ; T 2p ] ■ ■ ■ ',T 2 ffJ is an N y + 1-by-l vector. For 2 p th order accurate 
method, we have 

T kj = c kU (p+2 \x k ,yj,t), 0 < k < m — 1, 0 < j < N y , 

T 2p = 0, k> m 
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to the leading order, where m = 1,4, 6 , 8 ,11 for 2 p = 2,4, 6 , 8 ,10, Ck is a constant 
independent of h, t/( p+2 i := d ^ p +i and 0 is a zero vector of size N y + 1-by-l . 
We take the Fourier transform in y and Laplace transform in t , and obtain 


T% p = CkU( p+2 \xk,w,s), 0 < k < m — 1, 
T% p = 0, k > m. 


Similar to the one dimensional case, we only need to consider a vicinity of 
|s| = 0 with Re(s) > 0. By Lemma [foil we have Re(s+) > 0. The boundary 
part of Equation (1361) can be written in the matrix vector multiplication form 
C(s+,T)E UJ = hP+ 2 f^ B . 

We then need to investigate the invertibility of C(s+, t) and the dependence 
of \C~ 1 {s + ,t)T 2p \ on h. We note that Equation (fTTTTl) is analogous to the error 
equation m in one dimensional case, and the matrix-valued function C(-,r) 
is the same with only the argument changed. Therefore, based on the normal 
mode analysis for the one dimensional problem in t l 2. 3.111275751 for 2 p = 2,4,6 
and r > T 2 P the boundary error in the Fourier and Laplace space is bounded as 


\ 2 h < K u 


d p+2 

dx p+2 


U(0,v,s) 


2 

h 2p+ \ 


where K u depends on 77 but not h. Since the above estimate holds for every w, 
we sum them up and obtain 




d p+2 

dx p+2 


U(0,u>, 


s) 


2 

h 2p+ \ 


By using the Parseval’s relation in the Fourier space, we obtain 

2 


<K l 


p2iv 

d P+ ' 2 T -r, 

L 

dxP+ ^v>^‘) 


dy h 2p+ \ 


We then use Parseval’s relation in the Laplace space, and obtain 



dk >+ 2 
dx p+2 


U{0,y,t) 


2 

dydt. 


For 2 p = 2,4,6, the boundary error is bounded in 0(h p+2 ). We know from 
the energy estimate that the interior error is bounded in 0(h 2p ). Therefore, we 
conclude that for second, fourth and sixth order methods, the convergence rates 
are 2, 4 and 5 only if r > T 2 p. If r = T 2 p, the boundary system is singular. The 
convergence rate p + ^ given by the energy estimate is sharp. The numerical 
experiments in UHl agree with this conclusion. 

The above analysis can be carried out also for the Neumann problem and the 
interface problem (with conforming grid interfaces) in two dimensions. Com¬ 
paring with the corresponding one dimensional problem, the only difference in 
the error equation is that s is replaced by s+. By Lemma |5.II and the analysis 
thereafter, we conclude that the convergence result for the Neumann problem 
and the interface problem in two dimensions is the same as that for the corre¬ 
sponding one dimensional problem. For numerical experiments of these cases, 
we refer to m- 
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It is not straightforward to generalize the above analysis to a two dimensional 
problem with non-periodic boundary conditions at all boundaries, or to a two 
dimensional problem with a non-conforming grid interface. Many numerical 
experiments of these cases, however, show an agreement with the convergence 
rate that is conjectured from the analysis of the corresponding one dimensional 
problem. 


6 Numerical experiments 

In this section, we perform numerical experiments to verify the accuracy anal¬ 
ysis. We investigate how the accuracy of the numerical solution is affected by 
the large truncation error localized near the boundary and grid interface. In the 
analysis, we use the half line and half plane problem. However, in the numer¬ 
ical experiments, we use a bounded domain and impose boundary conditions 
on all boundaries weakly. We discretize in time using the classical fourth order 
Runge-Kutta method. The step size At in time is chosen so that the temporal 
error has a negligible impact of the spatial convergence result. The value of At 
is given in each numerical experiment. 

The L 2 error and maximum error are computed as the norm of the difference 
between the exact solution u ex and the numerical solution u h according to 

|| u h — u ex ||l 2 = \J h d (u h — u ex ) T (u h — u ex ), 

|| u h — u ex ||00 = max \u h — u ex |, 


where d is the dimension of the equation. The convergence rates are computed 

by 


q = log 





6.1 One dimensional wave equation 

We choose the analytical solution 


u = cos(107nr + 1) cos(107rf + 2), 0 < x < 1, 0 < t < 2, (37) 


and use it to impose the Dirichlet boundary conditions at x = 0 and x = 1. The 
analytical solution and its derivatives do not vanish at the boundaries at the 
final time t = 2. The L 2 error, the convergence rates in L 2 norm and maximum 
norm are shown in Table [4] The convergence rates behave as expected. With 
a 2p th (2p = 2,4,6) order accuracy method, the convergence rate is p + - in L 2 
norm if the penalty parameter r equals its limit. In this case, the pointwise 
error in the Laplace space is bounded in <D(h p ). We can hope for p th order 
convergence rate in the maximum norm, which is observed as well. With an 
increased penalty parameter, the accuracy of the numerical solution is improved 
significantly, and the convergence rate min(2 p,p + 2) is obtained. For second 
and fourth order accuracy cases, the accuracy difference between r = 1 . 2 t 2 P 
and r = 3t-2 P is almost invisible. For sixth order accuracy case, one can see the 
small improvement in accuracy between r = 1.2 t2 P and r = 3 t2 P . The step size 
in time is At = 0.1 h in all cases, except for the last two grid refinements of the 
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10" 7 3.97/3.98 
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Table 4: Convergence study for one dimensional wave equation with Dirichlet 


boundary conditions. 






Fourth Order 

Perturbed 


Sixth Order 

Perturbed 


N || u h - 

« ex ||L 2 <?l 2 

\\u h — U ex \\l, 2 

9l 2 

\\u h -u™\\ h2 

9l 2 \\u h -u ex \\i 

2 9L2 

51 1.84 

■ 10 -2 

1.46 ■ 10“ 2 


9.56 ■ 10“ 3 

9.56 ■ 10“ 3 


101 1.36 

■ 10" 3 3.76 

3.79 ■ 10" 3 

1.95 

1.94 ■ 10 -5 

8.94 7.40 ■ 10“ 5 

7.01 

201 8.17 

■ 10" 5 4.06 

6.62 ■ 10“ 4 

2.52 

3.04 ■ 10“ 7 

6.00 7.40 ■ 10 -6 

3.32 

401 4.88 

■ 10“ 6 4.07 

9.03 ■ 10“ 5 

2.88 

8.64 ■ 10“ 9 

5.14 5.12 ■ 10 -7 

3.85 

801 2.97 

. 10“ 7 4.04 

1.16 ■ 10 -5 

2.96 

1.20 • IQ” 10 

6.18 3.29 ■ 10 -8 

3.96 


Table 5: Convergence study for one dimensional wave equation with Neumann 
boundary condition. 


sixth order scheme with r = 3 t2 P where At = 0.05 h is used in order to observe 
the expected convergence behaviour. Note that these step sizes are much smaller 
than the step size restricted by the stability condition. With r = 1.2 t2 P , the 
numerical method is stable if A t/h < c, where the Courant number c = 1.3, 0.9 
and 0.6 for second, fourth and sixth order method. For an increased penalty 
parameter r = 3 t2 P , the corresponding Courant numbers are 0.8, 0.5 and 0.4. 

To test the Neumann problem, we use (1571) as the analytical solution and 
impose Neumann boundary condition at two boundaries x = 0 and x = 1. The 
L 2 error and the convergence rate in L 2 norm are shown in Table [5] in Column 
2 and 3 for the fourth order scheme. Fourth order convergence rate is clearly 
observed, which corresponds to the optimal accuracy order gain of two from 
the boundary truncation error. In the accuracy analysis in (J51 the optimal 
convergence rate relies on the fact that the truncation error vector is in the 
column space of C4jv(0). We add a dissipative term to the boundary block of 
the SBP operator so that the boundary truncation error vector is in the same 
magnitude but no longer in the column space of C4jv(0). The result is listed 
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Table 6: Convergence study for one dimensional wave equation with a grid 
interface. 


in Column 4 and 5 in Table [5] Now the gain in accuracy order is only one, 
and third order convergence rate is obtained as expected. Note that an energy 
estimate is also valid for the perturbed problem and guarantees stability. We 
repeat the same experiments with the sixth order scheme. The convergence 
rate is between five and six. We then add a dissipative term in the same way as 
for the fourth order scheme. The optimal accuracy order gain is lost, and the 
convergence rate drops to four. 

Next, we test how the accuracy and convergence are affected by the large 
truncation error localized near a grid interface. We choose the same analytical 
solution (El). The grid interface is located at x = 0.5 in the computational 
domain fi = [0,1], and the grid size ratio is 2 : 1. The step size in time is 
At = O.l/i, where h is the smaller grid size. Note that the analytical solution 
does not vanish at the grid interface. We use the SAT technique to impose the 
outer boundary condition weakly and choose the boundary penalty parameters 
strictly larger than its limit. The interface conditions are also imposed by the 
SAT technique. We use different interface penalty parameters to see how they 
affect the accuracy and convergence. The results are shown in Table 0 where 
N denotes the number of grid points in the coarse domain. 

According to the accuracy analysis in (JH for second and fourth order accu¬ 
racy methods with r = 72p, the expected convergence rates are 1.5 and 2.5 in 
L 2 norm. This is clearly observed in Table [6] We also observe in Table ED that 
the convergence rate in the maximum norm behaves as p. With an increased 
penalty parameter, a much better convergence result is obtained. For the second 
and fourth order methods, we get the optimal (second and fourth, respectively) 
order of convergence in both L 2 norm and maximum norm. For the sixth order 
scheme, we expect fifth order convergence in L 2 norm. The numerical exper- 
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101 
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Table 7: Convergence study for two dimensional wave equation. 


iments actually give a higher convergence result. This has also been observed 
for Schrodinger equation in mm- 

Remark The numerical experiments with 8 th and 10 t,t order SBP operators 
are also carried out. When the penalty parameter equals its lower limit, 4.5 and 
5.5 convergence rate in L 2 norm, and 4 and 5 convergence in maximum norm are 
clearly observed. With a larger penalty parameter, we expect the convergence to 
be 6 and 7, respectively. However, due to the high accuracy, the error decreases 
fast to the machine precision. Therefore, the expected convergence rates are 
not clearly observed. 

6.2 Two dimensional wave equation 

For two dimensional simulation, we choose the analytical solution 

u = cos(12x+l) cos(47ry+2) cos(\/l2 2 + (47r) 2 f+3), 0 < a; < 1, 0 < y < 1, 0 < t < 2, 

and use it to impose the Dirichlet boundary condition at x = 0 and x = 1. 
Periodic boundary condition is imposed at y = 0 and y = 1. The numbers of 
grid points are the same in x and y direction, and denoted by N in Table [7] 

The step size in time is At = 0.1 h. We clearly observe that with a 2 p th order 
method and r = T 2 p the convergence rates in L 2 norm are p+1 and in maximum 
norm are p. With an increased penalty parameter, the errors in L 2 norm are 
much smaller. The convergence rates are min(2 p,p + 2) in both L 2 norm and 
maximum norm. The results from the numerical experiments agree well with 
the analysis in £j5] 

7 Conclusion 

For the second order wave equation a stable numerical scheme does not automat¬ 
ically satisfy the determinant condition. We have considered stable SBP -SAT 
schemes for the Dirichlet, Neumann and interface problems. In all these cases 
the boundary truncation error is larger than the interior truncation error. We 
find that only the Dirichlet problem with penalty parameter greater than its 
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limit value satisfies the determinant condition, and in agreement with B QH, 
and a gain of two in accuracy order for the boundary error follows. Also in this 
case, the determinant condition is necessary for the optimal convergence rate. 

For all the other cases the determinant condition is not satisfied. For most 
of these cases there is a gain in accuracy of order one or even two, though 
the energy estimate only implies a gain of order 1/2. We show that a careful 
analysis of the effect of the boundary truncation error yields sharp estimates. 
We also show that the results from analysis in one space dimension are valid in 
two space dimensions. 

In particular we have found that to get the optimal gain in accuracy order 
for the Dirichlet problem and the grid interface problem, the penalty parameter 
should be chosen larger than its limit value given by stability analysis. However, 
we should not choose a very large penalty parameter since it does not decrease 
further the error in the solution but leads to a very small Courant number, thus 
an unnecessarily small step size in time. 


Appendix 1 


The matrix C+d+t) in £12.3.21 is 


C+(0, t ) 


'-3 

-1 

3 — 4\/3 
-1 

-122+48r 

ll 

59 

5 

2 

55 

85+12+3 

68 

59 

43 

43 

43 

43 

1 

37+3+3 

17 

0 

L 49 

49 

49 


The first four columns of the matrix C6 _d(0,t) in £ 12.3.31 are 


C6d(0,t) — 

s v ✓ 

column 1—4 


fc(r) 

1.823470407 

-4.136345752 

0.8614698016 

-0.0951377428 

-0.04002001476 


8.024525606 

2.35039818 

-4.137556867 

1.159078186 

-0.9156510516 

0.1875223549 


-8.215717879 

-1.867463026 

8.108447068 

-3.511693724 

1.987601033 

-0.3471267779 


2.979430728 

0.5704778157 

-4.615068241 

2.565973129 

-1.25102831 

0.1996244378 


where k(r) = 3.165067038r — 9.382128117, and the last two columns of the 
matrix Csd(0,t) in £12.3.31 are 


C6d(0,t) = 

Si v ✓ 

column 5—6 


3.368616929 + 0.02305363188i 
1.00246526 + 0.04694800264i 
-6.789853693- 0.2265966418i 
4.592688087+ 0.1965568935i 
-3.188300477- 0.2835125258i 
0.3280801113+ 0.1028755037i 


3.368616929- 0.02305363188i 
1.00246526 - 0.04694800264i 
-6.789853693 + 0.2265966418i 
4.592688087 - 0.1965568935i 
-3.188300477+ 0.2835125258i 
0.3280801113- 0.1028755037i 
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Appendix 2 

In M3.2.21 the boundary system for fourth and sixth order schemes are analyzed. 
We have for the fourth order scheme 



r 54 

59 

5 

11-4V3 


17 

17 

17 

17 

CW(0) = 

-1 

2 

-1 

-1 

4 

59 

55 

85+12+3 


43 

43 

43 

43 


1 

0 

1 

37+8+3 


L 49 

49 

49 


0 0 -0.0588 0 
0 0 0 0 
0 0 1.1860 0 ’ 
0 0 -0.0408 0 


and 


rp4:,Neu,B 

u 


For sixth order scheme, we have 


43 ' 
204 

w 

588 J 


UXXXX ( 0 ? 5 )• 


C6Jv(0) = 

column 1—4 


3.8056512077 

-1.0534129693 

0.6441780401 
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CW(0) = 
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1.0024652602 + 0.0469480026i 
-6.7898536930 - 0.2265966418i 
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-0.8104124471 + 0.0403172285i' 
1.0024652602- 0.0469480026i 
-6.7898536930 + 0.2265966418i 
4.5926880872- 0.1965568935i 
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0 0 0 -0.2598847290 0 0 

0 0 0 0.3269055745 0 0 

0 0 0 -1.7658674536 0 0 

0 0 0 1.8336101263 0 0 ’ 

0 0 0 -0.6935339173 0 0 
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rri6,NeU,B 
1 u 


'-0.3287481378' 

0.2200796359 

-0.5608447068 

0.2044006966 

-0.1710063053 

0.0514543047 


U x 


(0, s). 


Appendix 3 

In M.ll the Taylor expansion of C 2 i (5, r) at s = 0 is 

C 2 i(s,t) = C 21 (0, r) + sC' 2I (0, r) + 0(s 2 ), 
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where 


C 2 /( 0 , r) 
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In 94.21 the Taylor expansion of C 47 (s, r) at s = 0 is 
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